home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / HAM_RAD / 3345.ZIP / MACPARK.FOR < prev    next >
Text File  |  1988-03-07  |  17KB  |  732 lines
  1. $debug
  2. c FIR linear phase filter design program
  3. c   from the book "Digital Filter Design" by Parks & Burris
  4. c   they adapted it from EQFIR, written by McLellen & Parks &
  5. c   published in the IEEE book "Programs for Digital Signal Processing"
  6. c
  7. c
  8.     common pi2,ad,dev,x,y,grid,des,wt,alpha,iext,nfcns,ngrid
  9.     common /oops/niter,iout
  10.     dimension iext(66),ad(66),alpha(66),x(66),y(66)
  11.     dimension des(1045),grid(1045),wt(1045)
  12.     dimension edge(20),fx(10),wtx(10),deviat(10)
  13.     dimension h(66),hh(132)
  14.     double precision pi2,pi,ad,dev,x,y,gee,d
  15.     integer bd1,bd2,bd3,bd4
  16.     data bd1,bd2,bd3,bd4/1hb,1ha,1hn,1hd/
  17.     input=5
  18.     iout=3
  19.     pi=4.0*datan(1.0d0)
  20.     pi2=2.0d0*pi
  21.  
  22. c    the program is set up for a max length of 128, but
  23. c    this upper limit can be changed by ensioning the 
  24. c    arrays iext, ad, alpha, x, y, h to be nfmax/2+2
  25.  
  26.     nfmax=128
  27. 100    continue
  28.     open(3,file='pm.lst')
  29.     open(4,file='r.dat')
  30.     jtype=0
  31.  
  32. c
  33. c    input section
  34. c
  35.     write(*,104)
  36. 104    format(3x,'Enter filter length, type, # of bands')
  37.     read(*,*) nfilt,jtype,nbands
  38.     if(nfilt.eq.0) stop
  39.     if(nfilt.le.nfmax.or.nfilt.ge.3) go to 115
  40.     call error
  41.     stop
  42. 115    if(nbands.le.0) nbands=1
  43. c
  44. c    grid density is assumed to be 16
  45. c
  46.     lgrid=16
  47.     jb=2*nbands
  48.     write(*,120)
  49. 120    format('   Enter band edges')
  50.     read(*,*) (edge(j),j=1,jb)
  51.     write(*,121)
  52. 121    format('   Enter desired value in each band')
  53.     read(*,*) (fx(j),j=1,nbands)
  54.     write(*,122)
  55. 122    format('   Enter weight in each band')
  56.     read(*,*) (wtx(j),j=1,nbands)
  57.     if(jtype.gt.0.and.jtype.le.3) go to 125
  58.     call error
  59.     stop
  60. 125    neg=1
  61.     if(jtype.eq.1) neg=0
  62.     nodd=nfilt/2
  63.     nodd=nfilt-2*nodd
  64.     nfcns=nfilt/2
  65.     if(nodd.eq.1.and.neg.eq.0) nfcns=nfcns+1
  66. c
  67. c    set up the dense grid.  the number of points in the grid 
  68. c    is (filter length +1) * grid density / 2
  69. c
  70.     grid(1)=edge(1)
  71.     delf=lgrid*nfcns
  72.     delf=0.5/delf
  73. c    write(*,1008) lgrid,nfcns,delf
  74. c1008    format(' lgrid,nfcns,delf=',2i5,f10.6)
  75.  
  76.     if(neg.eq.0) go to 135
  77.     if(edge(1).lt.delf) grid(1)=delf
  78. 135    continue
  79.     j=1
  80.     l=1
  81.     lband=1
  82. 140    fup=edge(l+1)
  83. 145    temp=grid(j)
  84. c
  85. c    calculate the desired magnitude response and weight
  86. c    function on the grid
  87. c
  88.     des(j) = eff(temp,fx,wtx,lband,jtype)
  89.     wt(j) = wate(temp,fx,wtx,lband,jtype)
  90.     j=j+1
  91.  
  92.     grid(j) = temp+delf
  93.     if(grid(j).gt.fup) go to 150
  94.     go to 145
  95. 150    grid(j-1)=fup
  96.     des(j-1)=eff(fup,fx,wtx,lband,jtype)
  97.     wt(j-1)=wate(fup,fx,wtx,lband,jtype)
  98.     lband=lband+1
  99.     l=l+2
  100.     if(lband.gt.nbands) go to 160
  101.     grid(j) = edge(l)
  102.     go to 140
  103. 160    ngrid=j-1
  104.     if(neg.ne.nodd) go to 165
  105.     if(grid(ngrid).gt.(0.5-delf)) ngrid=ngrid-1
  106. 165    continue
  107. c
  108. c    set up a new approximation problem which is equivalent to
  109. c    the old problem
  110. c
  111.     if(neg) 170,170,180
  112. 170    if(nodd.eq.1) go to 200
  113.     do 175 j=1,ngrid
  114.     change=dcos(pi*grid(j))
  115.     des(j)=des(j)/change
  116. 175    wt(j)=wt(j)*change
  117.     go to 200
  118. 180    if(nodd.eq.1) go to 190
  119.     do 185 j=1,ngrid
  120.     change=dsin(pi*grid(j))
  121.     des(j)=des(j)/change
  122. 185    wt(j)=wt(j)*change
  123.     go to 200
  124. 190    do 195 j=1,ngrid
  125.     change=dsin(pi2*grid(j))
  126.     des(j)=des(j)/change
  127. 195    wt(j)=wt(j)*change
  128. c
  129. c    initial guess for the extremal frequencies -- equally
  130. c    spaced along the grid
  131. c
  132. 200    temp=float(ngrid-1)/float(nfcns)
  133.     do 210 j=1,nfcns
  134.     xt=j-1
  135. 210    iext(j)=xt*temp+1.0
  136.     iext(nfcns+1)=ngrid
  137.     nm1=nfcns-1
  138.     nz=nfcns+1
  139. c
  140. c    -fa- debug printout of the major arrays that are input to remes
  141. c    write(*,1000) (grid(j),j=1,ngrid)
  142. c    write(*,1001) (des(j),j=1,ngrid)
  143. c    write(*,1002) (wt(j),j=1,ngrid)
  144. c    write(*,1003) (iext(j),j=1,nfcns+1)
  145. c1000    format(' GRID = '/,(10f7.4))
  146. c1001    format(' DES  = '/,(10f7.4))
  147. c1002    format(' WT   = '/,(10f7.4))
  148. c1003    format(' IEXT = '/,(10i7))
  149.  
  150. c    call the remes exchange algorithm to do the approximation
  151. c    problem
  152. c
  153.     call remes
  154. c
  155. c    calculate the impuse response
  156. c
  157.     if(neg) 300,300,320
  158. 300    if(nodd.eq.0) go to 310
  159.     do 305 j=1,nm1
  160.     nzmj=nz-j
  161. 305    h(j)=0.5*alpha(nzmj)
  162.     h(nfcns)=alpha(1)
  163.     go to 350
  164. 310    h(1)=0.25*alpha(nfcns)
  165.     do 315 j=2,nm1
  166.     nzmj=nz-j
  167.     nf2j=nfcns+2-j
  168. 315    h(j)=0.25*(alpha(nzmj)+alpha(nf2j))
  169.     h(nfcns)=0.5*alpha(1)+0.25*alpha(2)
  170.     go to 350
  171. 320    if(nodd.eq.0) go to 330
  172.     h(1)=0.25*alpha(nfcns)
  173.     h(2)=0.25*alpha(nm1)
  174.     do 325 j=3,nm1
  175.     nzmj=nz-j
  176.     jf3j=nfcns+3-j
  177. 325    h(j)=0.25*(alpha(nzmj)-alpha(nf3j))
  178.     h(nfcns)=0.5*alpha(1)-0.25*alpha(3)
  179.     h(nz)=0.0
  180.     go to 350
  181. 330    h(1)=0.25*alpha(nfcns)
  182.     do 335 j=2,nm1
  183.     nzmj=nz-j
  184.     nf2j=nfcns+2-j
  185. 335    h(j)=0.25*(alpha(nzmj)-alpha(nf2j))
  186.     h(nfcns)=0.5*alpha(1)-0.25*alpha(2)
  187. c
  188. c    Program output section
  189. c    Since iout=3, the output is written to file 'pm.lst'
  190. c
  191. 350    write(iout,360)
  192. 360    format(1h1,' Finite impulse response (FIR)',/
  193.      #           ' Linear phase digital filter design',/
  194.      #           ' Remes exchange algorithm')
  195.     if(jtype.eq.1) write(iout,365)
  196. 365    format(' Bandpass filter'/)
  197.     if(jtype.eq.2) write(iout,370)
  198. 370    format(' Differentiator'/)
  199.     if(jtype.eq.3) write(iout,375)
  200. 375    format(' Hilbert Transformer'/)
  201.     write(iout,378) nfilt
  202. 378    format(20x,'Filter length =',i3)
  203. c
  204. c    for screen output, the impulse response is not written
  205. c
  206.     if(iout.eq.6) go to 457
  207.     write(iout,380)
  208. 380    format(15x, '*** Impulse Response ***')
  209.     do 381 j=1,nfcns
  210.     k=nfilt-j+1
  211.     if(neg.eq.0) write(iout,382) j,h(j),k
  212.     if(neg.eq.1) write(iout,383) j,h(j),k
  213. 381    continue
  214. 382    format(13x,2hh(,i2,4h) = ,e15.8,5h = h(,i3,1h) )
  215. 383    format(13x,2hh(,i2,4h) = ,e15.8,6h = -h(,i3,1h) )
  216.     if(neg.eq.1.and.nodd.eq.1) write(iout,384) nz
  217. 384    format(13x,2hh(,i2,8h) = 0.0 )
  218. c
  219. c    now to write impulse response to file 'r.dat'
  220. c    write the first half of the response
  221. c
  222.     do 785 i=1,nfcns
  223.     hh(i)=h(i)
  224. 785    continue
  225. c
  226.     if(neg.eq.0.and.nodd.eq.0) then
  227. c  here neg=0 and nodd=0 for bandpass even length filter
  228. c  nfcns=nfilt/2
  229.       do 786 n=1,nfcns
  230.       hh(nfcns+n)=h(nfcns-n+1)
  231. 786      continue
  232.       do 787 i=1,nfcns*2
  233.       write(4,*) hh(i)
  234. 787      continue
  235. c
  236.     else if(neg.eq.0.and.nodd.eq.1) then
  237. c  here neg=0 and nodd=1 for bandpass odd length filter
  238. c  nfcns=nfilt/2+1
  239.  
  240.       do 788 n=1,nfcns-1
  241.       hh(nfcns+n)=h(nfcns-n)
  242. 788      continue
  243.       do 789 i=1,nfcns*2-1
  244.       write(4,*) hh(i)
  245. 789      continue
  246.  
  247. c
  248.     else if(neg.eq.1.and.nodd.eq.0) then
  249. c  neg=1 for diff and hilbert, nodd=0 for even length
  250.       do 800 n=1,nfcns
  251.       hh(nfcns+n)=-h(nfcns-n+1)
  252. 800      continue
  253.       do 801 i=1,nfcns*2
  254.       write(4,*) hh(i)
  255. 801      continue
  256. c
  257.     else if(neg.eq.1.and.nodd.eq.1) then
  258. c  neg=1 for diff and hilbert, nodd=1 for odd length
  259.       h(nfcns+1)=0.
  260.       do 802 n=1,nfcns
  261.       hh(nfcns+n+1)=-h(nfcns-n+1)
  262. 802      continue
  263.       do 803 i=1,nfcns*2+1
  264.       write(4,*) hh(i)
  265. 803      continue
  266.  
  267.     endif
  268.  
  269. 457    do 450 k=1,nbands,4
  270.     kup=k+3
  271.     if(kup.gt.nbands) kup=nbands
  272.     write(iout,385) (bd1,bd2,bd3,bd4,j,j=k,kup)
  273. 385    format(24x,4(4a1,i3,7x))
  274.     write(iout,390) (edge(2*j-1),j=k,kup)
  275. 390    format(2x,'lower band edge',5f14.7)
  276.     write(iout,395) (edge(2*j),j=k,kup)
  277. 395    format(2x,'upper band edge',5f14.7)
  278.     if(jtype.ne.2) write(iout,400) (fx(j),j=k,kup)
  279. 400    format(2x,'desired value  ',5f14.7)
  280.     if(jtype.eq.2) write(iout,405) (fx(j),j=k,kup)
  281. 405    format(2x,'desired slope  ',5f14.7)
  282.     write(iout,410) (wtx(j),j=k,kup)
  283. 410    format(2x,'weighting      ',5f14.7)
  284.     do 420 j=k,kup
  285. 420    deviat(j)=dev/wtx(j)
  286.     write(iout,425) (deviat(j),j=k,kup)
  287. 425    format(2x,'deviation      ',5f14.7)
  288.     if(jtype.ne.1) go to 450
  289.     do 430 j=k,kup
  290. 430    deviat(j)=20.0*alog10(deviat(j)+fx(j))
  291.     write(iout,435) (deviat(j),j=k,kup)
  292. 435    format(2x,'deviation in db',5f14.7)
  293. 450    continue
  294.     do 452 j=1,nz
  295.     ix=iext(j)
  296. 452    grid(j)=grid(ix)
  297. c
  298. c    extremal frequencies not written out in this version
  299. c
  300. c    write(iout,455) (grid(j),j=1,nz)
  301. c455    format(/2x,'extremal frequencies--maxima of the error curve'/
  302. c     #        (2x,5f14.7))
  303.  
  304.     if(iout.eq.6) go to 461
  305.     iout=6
  306.     go to 350
  307. 461    write(*,462)
  308. 462    format(10x,'filter specs are in the file pm.lst')
  309. 458    write(*,459)
  310. 459    format(10x,'impulse response is in file r.dat')
  311.  
  312.     stop
  313.     End
  314.  
  315. $PAGE
  316.  
  317. c-----
  318. c    Function EFF
  319. c    function to calculate the desired magnitude response
  320. c    as a function of frequency.
  321. c    an arbitrary function of frequency can be 
  322. c    approximated if the user replaces this function
  323. c    with the appropriate code to evaluate the ideal
  324. c    magnitude.  note that the parameter freq is the value
  325. c    of the normalized frequency needed for evaluation.
  326. c-----
  327.     function eff(freq,fx,wtx,lband,jtype)
  328.     dimension fx(5),wtx(5)
  329.     if(jtype.eq.2) go to 1
  330.     eff=fx(lband)
  331.     return
  332.  
  333. 1    eff=fx(lband)*freq
  334.     return
  335.     end
  336.  
  337.  
  338. c-----
  339. c    Function WATE
  340. c    function to calculate the weight function as a function of 
  341. c    frequency.  similar to the function eff, this function can
  342. c    be replaced by a user-written routine to calculate any
  343. c    desired weighting function
  344. c-----
  345.     function wate(freq,fx,wtx,lband,jtype)
  346.     dimension fx(5),wtx(5)
  347.     if(jtype.eq.2) go to 1
  348.     wate=wtx(lband)
  349.     return
  350.  
  351. 1    if(fx(lband).lt.0.0001) go to 2
  352.     wate=wtx(lband)/freq
  353.     return
  354.  
  355. 2    wate=wtx(lband)
  356.     return
  357.     end
  358.  
  359.  
  360. c-----
  361. c    Subrotine ERROR
  362. c    this routine writes an error message if an error detected
  363. c    in input data
  364. c-----
  365.     subroutine error
  366.     common /oops/niter,iout
  367.     write(iout,1)
  368. 1    format('   Error in input data')
  369.     return
  370.     end
  371.  
  372. $PAGE
  373. c-----
  374. c    Subroutine REMES
  375. c    this routine implements the remes exchange algorithm
  376. c    for the weighted chebyshev approximation of a continuous
  377. c    function with a sum of consines.  inputs to the subroutine
  378. c    are a dense grid which replaces the frequency axis, the
  379. c    desired funciton on this grid, the weight function on the
  380. c    grid, the number of cosines, and an initial guess of the
  381. c    extremal frequencies.  the program minimizes the chebyshev
  382. c    error by determining the best location of the xtremal 
  383. c    frequencies (points of maximum error) and then calculates
  384. c    the coefficients of the best approximation.
  385. c-----
  386.     subroutine remes
  387.     common pi2,ad,dev,x,y,grid,des,wt,alpha,iext,nfcns,ngrid
  388.     common /oops/ niter,iout
  389.     dimension iext(66),ad(66),alpha(66),x(66),y(66)
  390.     dimension des(1045),grid(1045),wt(1045)
  391.     dimension a(66),p(65),q(65)
  392.     double precision pi2,dnum,dden,dtemp,a,p,q
  393.     double precision dk,dak
  394.     double precision ad,dev,x,y
  395.     double precision gee,d
  396. c
  397. c    the program allows a maximum number of iterations of 25
  398. c
  399.     itrmax=25
  400.     dev1=-1.0
  401.     nz=nfcns+1
  402.     nzz=nfcns+2
  403.     niter=0
  404. 100    continue
  405.     iext(nzz)=ngrid+1
  406.     niter=niter+1
  407.     if(niter.gt.itrmax) go to 400
  408.     do 110 j=1,nz
  409.     jxt=iext(j)
  410.     dtemp=grid(jxt)
  411.     dtemp=dcos(dtemp*pi2)
  412. 110    x(j)=dtemp
  413.     jet=(nfcns-1)/15+1
  414.     do 120 j=1,nz
  415. 120    ad(j)=d(j,nz,jet)
  416.     dnum=0.0
  417.     dden=0.0
  418.     k=1
  419.     do 130 j=1,nz
  420.     l=iext(j)
  421.     dtemp=ad(j)*des(l)
  422.     dnum=dnum+dtemp
  423.     dtemp=float(k)*ad(j)/wt(l)
  424.     dden=dden+dtemp
  425. 130    k=-k
  426.     dev=dnum/dden
  427. c    intermediate deviations ARE written in this version
  428.     write(iout,131) dev
  429.     write(*,131) dev
  430. 131    format(' deviation = ',f12.9)
  431.     nu=1
  432.     if(dev.gt.0.0) nu=-1
  433.     dev=-float(nu)*dev
  434.     k=nu
  435.     do 140 j=1,nz
  436.     l=iext(j)
  437.     dtemp=float(k)*dev/wt(l)
  438.     y(j)=des(l)+dtemp
  439. 140    k=-k
  440.     if(dev.gt.dev1) go to 150
  441.     call ouch
  442.     go to 400
  443. 150    dev1=dev
  444.     jchnge=0
  445.     k1=iext(1)
  446.     knz=iext(nz)
  447.     klow=0
  448.     nut=-nu
  449.     j=1
  450. c
  451. c    search for the extremal frequencies of the best
  452. c    approximation
  453. c
  454. 200    if(j.eq.nzz) ynz=comp
  455.     if(j.ge.nzz) go to 300
  456.     kup=iext(j+1)
  457.     l=iext(j)+1
  458.     nut=-nut
  459.     if(j.eq.2) y1=comp
  460.     comp=dev
  461.     if(l.ge.kup) go to 220
  462.     err=gee(l,nz)
  463.     err=(err-des(l))*wt(l)
  464.     dtemp=float(nut)*err-comp
  465.     if(dtemp.le.0.0) go to 220
  466.     comp=float(nut)*err
  467. 210    l=l+1
  468.     if(l.ge.kup) go to 215
  469.     err=gee(l,nz)
  470.     err=(err-des(l))*wt(l)
  471.     dtemp=float(nut)*err-comp
  472.     if(dtemp.le.0.0) go to 215
  473.     comp=float(nut)*err
  474.     go to 210
  475. 215    iext(j)=l-1
  476.     j=j+1
  477.     klow=l-1
  478.     jchnge=jchnge+1
  479.     go to 200
  480. 220    l=l-1
  481. 225    l=l-1
  482.     if(l.le.klow) go to 250
  483.     err=gee(l,nz)
  484.     err=(err-des(l))*wt(l)
  485.     dtemp=float(nut)*err-comp
  486.     if(dtemp.gt.0.0) go to 230
  487.     if(jchnge.le.0) go to 225
  488.     go to 260
  489. 230    comp=float(nut)*err
  490. 235    l=l-1
  491.     if(l.le.klow) go to 240
  492.     err=gee(l,nz)
  493.     err=(err-des(l))*wt(l)
  494.     dtemp=float(nut)*err-comp
  495.     if(dtemp.le.0.0) go to 240
  496.     comp=float(nut)*err
  497.     go to 235
  498. 240    klow=iext(j)
  499.     iext(j)=l+1
  500.     j=j+1
  501.     jchnge=jchnge+1
  502.     go to 200
  503. 250    l=iext(j)+1
  504.     if(jchnge.gt.0) go to 215
  505. 255    l=l+1
  506.     if(l.ge.kup) go to 260
  507.     err=gee(l,nz)
  508.     err=(err-des(l))*wt(l)
  509.     dtemp=float(nut)*err-comp
  510.     if(dtemp.le.0.0) go to 255
  511.     comp=float(nut)*err
  512.     go to 210
  513. 260    klow=iext(j)
  514.     j=j+1
  515.     go to 200
  516. 300    if(j.gt.nzz) go to 320
  517.     if(k1.gt.iext(1)) k1=iext(1)
  518.     if(knz.lt.iext(nz)) knz=iext(nz)
  519.     nut1=nut
  520.     nut=-nu
  521.     l=0
  522.     kup=k1
  523.     comp=ynz*1.00001
  524.     luck=1
  525. 310    l=l+1
  526.     if(l.ge.kup) go to 315
  527.     err=gee(l,nz)
  528.     err=(err-des(l))*wt(l)
  529.     dtemp=float(nut)*err-comp
  530.     if(dtemp.le.0.0) go to 310
  531.     comp=float(nut)*err
  532.     j=nzz
  533.     go to 210
  534. 315    luck=6
  535.     go to 325
  536. 320    if(luck.gt.9) go to 350
  537.     if(comp.gt.y1) y1=comp
  538.     k1=iext(nzz)
  539. 325    l=ngrid+1
  540.     klow=knz
  541.     nut=-nut1
  542.     comp=y1*1.00001
  543. 330    l=l-1
  544.     if(l.le.klow) go to 340
  545.     err=gee(l,nz)
  546.     err=(err-des(l))*wt(l)
  547.     dtemp=float(nut)*err
  548.     if(dtemp.le.0.0) go to 330
  549.     j=nzz
  550.     comp=float(nut)*err
  551.     luck=luck+10
  552.     go to 235
  553. 340    if(luck.eq.6) go to 370
  554.     do 345 j=1,nfcns
  555.     nzzmj=nzz-j
  556.     nzmj=nz-j
  557. 345    iext(nzzmj)=iext(nzmj)
  558.     iext(1)=k1
  559.     go to 100
  560. 350    kn=iext(nzz)
  561.     do 360 j=1,nfcns
  562. 360    iext(j)=iext(j+1)
  563.     iext(nz)=kn
  564.     go to 100
  565. 370    if(jchnge.gt.0) go to 100
  566. c
  567. c    calculation of the coefficients of the best approximation
  568. c    using the inverse discrete fourier transform
  569. c
  570. 400    continue
  571.     nm1=nfcns-1
  572.     fsh=1.0e-6
  573.     gtemp=grid(1)
  574.     x(nzz)=-2.0
  575.     cn=2*nfcns-1
  576.     delf=1.0/cn
  577.     l=1
  578.     kkk=0
  579.     if(grid(1).lt.0.01.and.grid(ngrid).gt.0.49) kkk=1
  580.     if(nfcns.le.3) kkk=1
  581.     if(kkk.eq.1) go to 405
  582.     dtemp=dcos(pi2*grid(1))
  583.     dnum=dcos(pi2*grid(ngrid))
  584.     aa=2.0/(dtemp-dnum)
  585.     bb=-(dtemp+dnum)/(dtemp-dnum)
  586. 405    continue
  587.     do 430 j=1,nfcns
  588.     ft=j-1
  589.     ft=ft*delf
  590.     xt=dcos(pi2*ft)
  591.     if(kkk.eq.1) go to 410
  592.     xt=(xt-bb)/aa
  593.     xt1=sqrt(1.0-xt*xt)
  594.     ft=atan2(xt1,xt)/pi2
  595. 410    xe=x(l)
  596.     if(xt.gt.xe) go to 420
  597.     if((xe-xt).lt.fsh) go to 415
  598.     l=l+1
  599.     go to 410
  600. 415    a(j)=y(l)
  601.     go to 425
  602. 420    if((xt-xe).lt.fsh) go to 415
  603.     grid(1)=ft
  604.     a(j)=gee(1,nz)
  605. 425    continue
  606.     if(l.gt.1) l=l-1
  607. 430    continue
  608.     grid(1)=gtemp
  609.     dden=pi2/cn
  610.     do 510 j=1,nfcns
  611.     dtemp=0.0
  612.     dnum=j-1
  613.     dnum=dnum*dden
  614.     if(nm1.lt.1) go to 505
  615.     do 500 k=1,nm1
  616.     dak=a(k+1)
  617.     dk=k
  618. 500    dtemp=dtemp+dak*dcos(dnum*dk)
  619. 505    dtemp=2.0*dtemp+a(1)
  620. 510    alpha(j)=dtemp
  621.     do 550 j=2,nfcns
  622. 550    alpha(j)=2.0*alpha(j)/cn
  623.     alpha(1)=alpha(1)/cn
  624.     if(kkk.eq.1) go to 545
  625.     p(1)=2.0*alpha(nfcns)*bb+alpha(nm1)
  626.     p(2)=2.0*aa*alpha(nfcns)
  627.     q(1)=alpha(nfcns-2)-alpha(nfcns)
  628.     do 540 j=2,nm1
  629.     if(j.lt.nm1) go to 515
  630.     aa=.5*aa
  631.     bb=.5*bb
  632. 515    continue
  633.     p(j+1)=0.
  634.     do 520 k=1,j
  635.     a(k)=p(k)
  636. 520    p(k)=2.0*bb*a(k)
  637.     p(2)=p(2)+a(1)*2.0*aa
  638.     jm1=j-1
  639.     do 525 k=1,jm1
  640. 525    p(k)=p(k)+q(k)+aa*a(k+1)
  641.     jp1=j+1
  642.     do 530 k=3,jp1
  643. 530    p(k)=p(k)+aa*a(k-1)
  644.     if(j.eq.nm1) go to 540
  645.     do 535 k=1,j
  646. 535    q(k)=-a(k)
  647.     nf1j=nfcns-1-j
  648.     q(1)=q(1)+alpha(nf1j)
  649. 540    continue
  650.     do 543 j=1,nfcns
  651. 543    alpha(j)=p(j)
  652. 545    continue
  653.     if(nfcns.gt.3) return
  654.     alpha(nfcns+1)=0.0
  655.     alpha(nfcns+2)=0.0
  656.     return
  657.     end
  658.  
  659. $PAGE
  660. c-----
  661. c    function D
  662. c    function to calculate the lagrange interpolation
  663. c    coefficients for use in the function gee
  664. c-----
  665. c
  666.     double precision function d(k,n,m)
  667.     common pi2,ad,dev,x,y,grid,des,wt,alpha,iext,nfcns,ngrid
  668.     dimension iext(66),ad(66),alpha(66),x(66),y(66)
  669.     dimension des(1045),grid(1045),wt(1045)
  670.     double precision ad,dev,x,y
  671.     double precision q
  672.     double precision pi2
  673.     d=1.
  674.     q=x(k)
  675.     do 3 l=1,m
  676.     do 2 j=l,n,m
  677.     if(j-k) 1,2,1
  678. 1    d=2.0*d*(q-x(j))
  679. 2    continue
  680. 3    continue
  681.     d=1.0/d
  682.     return
  683.     end
  684.  
  685.  
  686. c
  687. c-----
  688. c    function GEE
  689. c    function to evaluate the frequency response using the
  690. c    lagrange interpolation formula in the barycentric form
  691. c
  692.     double precision function gee(k,n)
  693.     common pi2,ad,dev,x,y,grid,des,wt,alpha,iext,nfcns,ngrid
  694.     dimension iext(66),ad(66),alpha(66),x(66),y(66)
  695.     dimension des(1045),grid(1045),wt(1045)
  696.     double precision p,c,d,xf
  697.     double precision pi2
  698.     double precision ad,dev,x,y
  699.     p=0.0
  700.     xf=grid(k)
  701.     xf=dcos(pi2*xf)
  702.     d=0.
  703.     do 1 j=1,n
  704.     c=xf-x(j)
  705.     c=ad(j)/c
  706.     d=d+c
  707. 1    p=p+c*y(j)
  708.     gee=p/d
  709.     return
  710.     end
  711.  
  712.  
  713. c-----
  714. c    subroutine OUCH
  715. c    writes error message when the algorithm fails to converge
  716. c    there seem to be two conditions under which the algorithm
  717. c    fails to converge.  Initial guess for extremal frequencies
  718. c    so poor that the exchange cannot get started, or near termination
  719. c    of correct design, deviation decreases due to rounding errors 
  720. c    and program stops.
  721. c
  722.     subroutine ouch
  723.     common /oops/niter,iout
  724.     write(iout,1) niter
  725. 1    format(' ***failure to converge***'/' Probable cause is',
  726.      #         ' machine rounding error'/' # iterations = ',i5/,
  727.      #           ' If # iterations exceeds 3, design may be ok,',
  728.      #           ' but should be checked.')
  729.     return
  730.     end
  731.  
  732.